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INTRODUCTION 

The problem of passive scalar transport by turbulent 
flows has received much attention lately. The progress 
achieved has been made possible mainly by the introduc- 
tion of the Kraichnan model M. Within the model the 
turbulent velocity statistics is believed to be Gaussian, 
scale invariant in space and S correlated in time, which 
allows one to write closed equations on the correlation 
functions of the scalar. Such a velocity has only a few 
rough features in common with realistic flows, which are 
intermittent and have a finite scale-dependent correlation 
time, contrary to what is assumed in the model. Never- 
theless, it seems that many interesting properties of the 
statistics are inherent in the dynamics, rather than due 
to the intermittency of the velocity statistics itself. Un- 
usually for the turbulence theory numerous results have 
been obtained analytically using the Kraichnan model. 

Having reached an understanding of this model, it is 
then natural to generalize its results, passing to more 
realistic flows. However, due to the complicated inter- 
play between spatial and temporal properties of the ve- 
locity, one encounters various difficulties in introducing 
a meaningful velocity field with a finite correlation time. 
The only case where this was easily done is the so called 
Batchelor regime ||, where the spatial structure of the 
velocity is rather simple, and therefore one can separate 
space and time dependencies. It appears in the limit of 
large Prandtl numbers, which is the ratio of the fluid vis- 
cosity to the diffusivity of the transported quantity. In 
studying advection below the viscous length, the corre- 
lation functions of the velocity are smooth functions of 
space, which allows one to introduce an effective descrip- 
tion with v a — a afj {t)rp gj. In this way time and space 
become completely separated. 

The Batchelor limit is well studied if a has a zero corre- 
lation time and its statistics is Gaussian [p|-p^[; that is, in 
the framework of the Kraichnan model. Certain results 
have been derived for arbitrary statistics of a |l]-||,[l3| . 

Our aim here is to investigate the degree of universal- 
ity of the passive scalar statistics for arbitrary temporal 
correlations of the velocity. We utilize the close rela- 



tion between the statistics of Lagrangian trajectories in 
a turbulent flow and the statistics of the passive scalar. 
Therefore, it seems reasonable to separate the problems, 
first investigating the Lagrangian motion and next apply- 
ing the results to particular problems. For the Batchelor 
regime only a few degrees of freedom characterize the La- 
grangian dynamics, which makes the problem solvable. 

The plan of this paper is as follows. First we pass to 
the comoving reference frame in the equation for a passive 
scalar, which allows us to consider the Lagrangian map- 
ping as an affine transformation, characterized by a ran- 
dom matrix. After its probability distribution function 
(PDF) is found, we consider several particular examples 
of the scalar statistics both for the decaying and forced 
turbulence. We show that the statistics of the scalar can 
be found by integration of the distribution function with 
a kernel, depending on the problem in question. 

I. GENERAL RELATIONS 

Advection of a passive scalar $ by incompressible ve- 
locity field v is described by the equation 

d t ti + 0,V)tf - kV 2 i? = 0, (1.1) 

where k is the molecular diffusivity. We shall be inter- 
ested in the limit of small but finite k. In the case of con- 
tinuous injection of the scalar, one sh ould add a source 
<t>{t,r) into the right-hand side of Eq. fll.l| ). 

Let us consider a blob of the scalar having a size L 
much smaller than the viscous length of the velocity. The 
variation of the velocity on the scale of the blob is much 
smaller than the large homogeneous velocity transferring 
the blob as a whole. To account for a slow variation of the 
form of the blob due to relative motion of the particles, 
it is natural to pass to the reference frame moving with 
the velocity of a particle within the blob |IJ,|l5|. Since 
the velocity is a smooth function on the scale of the blob, 
it can be expanded in a Taylor series, thus leading to the 
equation 

d t & + o- a f3r p V a $ - nV 2 d = . (1.2) 
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Here o a p (t) is the matrix of the velocity derivatives taken 
at the chosen Lagrangian point. Incompressibility im- 
plies a aa = 0. For turbulent flows a should be regarded 
as a random matrix, having a finite correlation time r, 
which is the Lagrangian correlation time of the velocity. 

The complete information about the Lagrangian flow, 
defined by dtR = <tR, is contained in the matrix W, 
satisfying 

d t W = aW , W(0) = 1 . 

This generates an affine transformation of space points, 
so that a vector R transforms as R(t) = W(t)R(0). The 
volume conservation guarantees det W — 1. The motion 
of the particles of the scalar differs from that of the space 
points due to the nonzero diffusivity. To investigate this 
motion wc introduce the " inertia tensor" of the blob [|l4| , 



i a p 



1 

2N 



drr a rj3-d{t,r) 



(1.3) 



where N = jdri3(t,r) is the number of particles of the 
scalar. It is e asy to check that N is conserved by the 
full equation ( |l.2| ). It turns out that / contains all the 
necessary information, and will appear in the following 
sections as the result of formal calculations. The tensor 
I satisfies the closed dynamical equation 



d t I = k + al + la 1 



(1.4) 



The initial condition depends on the form of the ini- 
tial blob, generally I a p ~ L 2 . We shall see that for 
problems with spatial isotropy it is enough to consider 
lap(0) — L 2 5 a 0. One can check that / can be expressed 
via W in a way that is nonlocal in time. Since / is sym- 
metric and incorporates the diffusion, instead of working 
with W it will be more convenient for us to work directly 
with Eq. (jul) . 

The dynamics of the symmetric matrix / can be sepa- 
rated into the nontrivial essential dynamics of its eigen- 
values and the trivial dynamics of the angular degrees of 
freedom. It is thus natural to reformulate the dynamics 
for the eigenvalues directly, excluding irrelevant angular 
degrees of freedom. For the case that is S correlated in 
time, this can be done exactly, resulting in the Calogero- 
Sutherland model ||. We shall show that for a finite 
correlation time of a the angular degrees of freedom can 
also be effectively excluded. The reason for this is that 
only the large-time dynamics of the eigenvalues is impor- 
tant for our purposes so that in many respects (but not 
all) the matrix a appears to be 5 correlated in time. 

Before wc proceed with the derivation, it is useful 
to understand qualitatively the typical dynamics of a 
blob. If the amplitude of the velocity fluctuations is large 
enough (the precise condition will be formulated below), 
the term k on the right-hand side of Eq. (1.4) can be dis- 
regarded during the initial stage of the evolution. Then 
one can make sure that I coincides with L 2 WW T . Ac- 
cording to the Oseledets theorem [161, at large enough 



times the logarithms of the eigenvalues of the latter ma- 
trix are asymptotically equal to 2Xit. The Lyapunov ex- 
ponents Ai, . . ., X d do not depend on a particular real- 
ization of ct; hence they are important characteristics of 
the system. 

The above implies that the directions corresponding 
to positive and negative Lyapunov exponents will grow 
or decrease correspondingly, and that the blob will be- 
come an ellipsoid with the length of its main axis chang- 
ing as exp(A,i). The orientation of the ellipsoid can be 
arbitrary. The smallest dimension will decrease exponen- 
tially with the rate Ad, until at t |Ad| _1 In (L 2 \Xd\/K) 
it reaches a scale r^a = y/n/\Xd\, where the diffusive 
spreading of particles makes further contraction impos- 
sible. Later, the smallest direction will fluctuate around 
rdif. This will not affect other directions, that will con- 
tinue to change according to their Lyapunov exponents. 
In order to have a wide separation of the scales L and rdif , 
one should require a large value of the Peclet number 



Pe = L/r dif = Ly/\\ d \/n. 



(1.5) 



This ensures that the time needed to reach the diffusion 
scale is large, so that the above arguments are valid. 

Apart from the typical event described here, we shall 
also need the distribution of all outcomes. This is the 
aim of Sec. ||. Although it is not difficult to work with 
the arbitrary dimensionality of space, we shall consider 
the physical dimensionalities d — 2 and 3 only. 



II. STATISTICS OF / 

To separate the angular degrees of freedom from the 
radial ones, it is natural to represent / as follows: 



I = R T AR . 



(2.1) 



Here R is an orthogonal matrix composed of the eigen- 
vectors of /, and A is a diagonal matrix with the eigen- 
values e 2pi , . . . e 2/>d along the diagonal (we believe that 
the eigenv alue s are ordered, so that p\ > p2 > • •• > Pd)- 
Equation (1.4) becomes 



dtPi = d~ii + - exp(-2pi), a = RaR T , (2.2) 



d t R = nR , a, = 



g 2pi _ g2pj 



(2.3) 



We do not as sum e a s umm ation over the repeating in- 
dices in Eqs. (2_^2) and (^3). This system of equations is 
not very useful for analyzing the general case. However, 
one can notice that if during the evolution the eigenvalues 
become widely separated, that is pi ^> • • • 3> pa, the sys- 
tem is greatly simplified. In this case the antisymmetric 
matrix f2 becomes 



a- 



i < k 
i > k . 



(2.4) 
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and due to Eqs. (2.2.2.3) the dynamics of the angu- 
lar degrees of fr eedo m is independent on the eigenvalues. 
Therefore, Eq. (2.2) can be resolved: 

Pi =Poi+ I dt'd u {t') (2.5) 
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1 + ne- 2pm J dt' exp < 


-2jT dt'a u (t") 


\] 









■In 



Here poi are some constants of the order of unity that 
should be determined by matching with the initial pe- 
riod of separation of the eigenvalues. 

Integrals in Eq. ( |2.5| ) determine the dynamics of pj. 
We consider times much larger than f , the correlation 
time of a, which is generally less than or of the order 
of r. The form of the probability distribution function 
of p is different if we consider f <C t < | A^ | ~~ 1 In Pe or 
t ^ | Arf | 1 InPe, depending on whether the diffusion has 
started to be relevant or not. If the former case is c onsid- 
ered, one can disregard the second term in Eq. (|2.5|), thus 
obtaining pi w f Q dt' du(t'). We recognize the case of the 
central limit theorem. However, the Gaussian distribu- 
tion describes only the bulk of the most probable events 
leaving rare events out of the domain of its validity. We 
shall need a more general expression jL7| which can be 
derived from the following considerations. The integrals 
can be considered as sums of a large number n w t/f 
of independent identically distributed random variables. 
Thus we investigate the distribution of X given by 



X 



(2.6) 



V oc exp 

V oc exp 



-tS 2 



-ts* 



Pi - 
t 



pi - Xit p 2 — X 2 t 



t 



t 



Without loss of the generality we can assume that 
(xi) = 0. If the generating function (ex.p[iyxi\) of each 
x is exp[— s(y)], then X has the generating function 
exp[— ns(y)]. To find the distribution function of X, one 
should make the inverse Fourier transform 



P(X) 



dy 
2tt 



exp[iyA - ns(y)] 



At large n this integral can be calculated in the saddle- 
point approximation. Writing the extremum condition, 
we see that y e xtr is a function of the argument X/n, 
which implies V(X) cx exp[— nS(X/n)]. For X <C n 
one can expand 5 in the Taylor series and obtain V oc 
exp[— X 2 /(2nA)}. Here A is the variation of Xi. This 
is nothing but the central limit theorem. On the other 
hand, if we increase n, keeping the ratio X/n a constant 
of the order of unity, we can assert that InV oc — n. This 
has a simple interpretation. If X is of the order of n, 
only realizations where most of of the same sign 

contribute. Therefore we can model the situation by the 
binomial random process, which gives just the above re- 
sult. 

If we replace sum ( |2.6P by the integral J dt'x(t') of a 
random function x over time t much larger than the cor- 
relation time, we should only note that the characteristic 
function of X is proportional to exp[— ts(y)], and then 
proceed as above. We used the fact that the character- 
istic function is an exponent of the cumulant generating 
function. The derivation is easily generalized for several 
quantities. Thus the distribution functions in d = 2 and 
3 are given by the formulas 



Q(Pi)8(Pi +P2) ■ 



0{pi - P2)0(p2 - p?,)8(pi + P2 + P3) ■ 



(2.7) 
(2.8) 



Here 52(21) and S^(x\,X2) are some functions depend- 
ing on the details of the statistics of a. In the 8 corre- 
lated case one can find the explicit expression for 52,3 
H[l3| (see also Appendix |a|) . The constants Ai are noth- 
ing but the Lyapunov exponents, which are expressed via 
the statistics of a in the following way (cf . Ref . p6| ) 



(2.9) 



To have a self-consistent picture, we should assume that 
the spectrum of the Lyapunov exponents is nondegener- 
ate, that is Ai > Aj+j.. Physically nondegeneracy of the 
Lyapunov exponents means that a blob is unstable with 
respect to the fluctuations, leading to a separation of the 
lengths of its sides. Noticing that the Lagrangian point 
r = is a saddle point for the incompressible flow, one 
can easily verify that the strain directions corresponding 
to the further elongation of the blob prevail. Therefore, 
during a time r of approximately constant strain, the 



blob will be on average further elongated. 

At t ^> (Xi — Ai + i) _1 we can disregard effects orig inat- 
i ng from the boundary pi — Pi+i. Equations (2/7) and 
( p.8[ ) are not valid in a narrow region near the boundary 
which has a width of the order of unity. Since it is much 
smaller than Xit, we can use the step function 9 to model 
the form of the PDF near the boundary. 

Due to the incompressibility condition the exponents 
satisfy 5Z i=1 Aj = 0. Then, in order to have a spectrum 
that is nondegenerate in d = 2, one should only require 
that Ai > 0, which is the same saying that trajectories 
diverge exponentially. In d — 3, it is necessary to supply 
some information about A2 . If the statistics of a is sym- 
metric with respect to time reversion, then A2 = [p_3| - 
However, if this is not the case, it is generally nonzero. 
In Appendix [b| we find the expression for A2 if the cor- 
relation time of a is small, which shows that its sign is 
generally arbitrary. 
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The form of the functions 52,3 depends on particu- 
lar details of the statistics of a. However, it is possible 
to make two general statements about these functions. 
First, one can assert that at small x the expansion 



2C_ 



11 



C 22 x\ - 2C 12 x x x 2 + Cnx\ 
%C\\C 22 — C 12 ) 



is valid, reproducing the central limit theorem. The con- 



stants Cij are defined as 



dt' 



Here ((...)) stands for irreducible correlation function. 
The integrals should be calculated over an interval, much 
larger then the correlation time of a. Note that the con- 
dition of incompressibility ensures that X^jLi &ij = 0- 

When Xi are of the order of unity, the functions 
52,3 have no singularities and change smoothly. The 
quadratic expansion of S 2t 3 is valid as long as 



\pi - A 4 t| < t/r, 



(2.10) 



m d = 3. 
1 lnPc. The diffusion is irrcl- 



where f is the correlation time of a. In the 5-correlated 
case it holds everywhere (Appendix ^J). 

The normalization of V is determined by the quadratic 
part of 5*2,3, since most of the probability is concen- 
trated at | pi — Xit\ ~ yj dj t <C t. One can find 
the normalization factors (2nCiit)~ 1 ^ 2 in d = 2 and 

[47r 2 t 2 (C n C 22 -Cf 2 )]- 1/2 : 

Now consider t > | 
evant for pi having a non- negative Lyapunov exponent. 
However, there is a finite probability, increasing with t, 
that pd reaches the diffusion scale. This requires an ac- 
count of the last term on the right-hand side of Eq. ( |2.2| ) 
or Eq. (2.5). The diffusion will not allow pd to decrease 
much below ln(|Ad|/«0- On the contrary, negative Ad will 
prevent it from increasing. As a result, the correspond- 
ing pi will be distributed stationarily around the value 
ln(«;/|Ad|). Relaxation times associated with this distri- 
bution are diffusion independent and thus are much less 
then t. On the other hand, p having non-negative Lya- 
punov exponents are the integrals over the whole evolu- 
tion time i, so that their values at time t are not sensitive 
to the last period of evolution with duration of the or- 
der of the relaxation time of pd- This means that fixing 
their values at time t ^S> |Ad| _1 mPe will not affect the 
distribution of pd, and the whole probability distribution 
function V is factorized (cf. Refs. (HH). In d = 2 we can 
write 



V oc exp < —t S; 



Pi - Ait 
t 



V s t{p 2 



(2.11) 



Here V s t is the stationary distribution of p 2 . In d = 3 the 
situation is more complicated. While A3 is always nega- 
tive, A2 can be both positive and negative. The form of 
the PDF will be different for these two cases. If A2 > 0, 



V oc exp < —t S; 



pi - Ait p 2 - A 2 t 



t 



t 



V st (p 3 ). (2.12) 



Since pd is independent of t he r est of as the functions 
62,3 are the same as in Eqs. ( [2.7[ ) and (2.8). 

If A2 is negative, then at t ^ | A2 1 1 ln(| A2 1/«) the 
distribution over p 2 will also become steady and con- 
centrated near ln(| A2 which by order of magnitude 
is equal to ln(| A3 Therefore, our assumption that 
Pi 3> P3 is incorrect. Still pi 3> P2,3, and the equation 
for pi is separated from the other variables. Then the 
distribution function is equal to 



V oc exp < —t S; 



Pi - Ait 
t 



V st {p 2 ,pn) 



(2.13) 



where S3 is related to S3 by exp(— S3) oc J dp 2 exp(— S3). 

Finally, let us note that since the configuration space 
SO(d) of the rotation matrix R [sec Eq. (jOj] is finite, 
and since there is no preferred direction in space, at large 
t the matrix is distributed uniformly over the sphere. 

The basic result obtained above is the special scaling 
form of the probability density functions. It is this uni- 
versal form which lies in the origin of the results derived 
below. 



III. DECAYING TURBULENCE 



As a first application, let us consider decay of a passive 
scalar d. The problem is posed as follows: given a ran- 
dom distribution of the scalar density $0 at t = 0, find 
its statistics at t > 0. In the framework of the Kraich- 
nan model the single-point statistics was considered by 
Son H , who obtained the following long-time asymptotic 
behavior: 



(m,0)\ a ) ocexp(- 7a t). 



(3.1) 



where j a in d = 3 is equal to a(6 — a)Z?/4 for < a < 3 
and 9D /4 otherwise (D is a parameter characterizing the 
strength of the fluctuations of a). The same decay law 
has been claimed for the gradients of the scalar. Here 
we consider the problem for arbitrary correlation time 
of a. Our consideration shows that due to the above- 
mentioned special form of Eqs. (|2.1l| - |2~T3| ) the law (3.1) 
is valid for an arbitrary statistics of a both for the single- 
point value of the scalar and its gradient. In the 5- 
correlated limit we obtain a result for 7 Q different from 
that of Ref. ||. The results also show that the basic 
assumption of Ref. ]18| is incorrect. 

The following qualitative picture, supported by the 
calculations presented below, explains the decay. First, 
consider a single blob initially having a characteristic 
size L and containing N particles of the scalar. As ve- 
locity stretches the blob, the number of particles does 
not change, contrary to the volume of the blob. At 
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t ^ |Ad|~ 1 ln(|Ad|^ 2 /«) the dimensions of the blob with 
negative Lyapunov exponents are frozen at rdif , while the 
rest keep growing exponentially, resulting in an exponen- 
tial growth of the total volume of the blob. Since the 
volume is proportional to v det / = exp(^p;), one has 
oc (exp (— a^2 pi)), where the averaging should be 
done with the help of the PDF discussed above, that 
is (|$| Q ) oc J dpV(p) exp (—aJ2 Pi)- The result is de- 
termined by a compromise between two competing fac- 
tors. While the averaged quantity exp(— aJ2 Pi) favors 
smaller values of ^2pi, the maximum of the probability 
is attained when each growing pi is equal to \t. Obvi- 
ously, for larger a the volume acquires more importance, 
so that the main contribution is made by smaller blobs, 
which are less probable but have larger concentration of 
the scalar. So, at small a, the deviation from the av- 
erage growth Xt is small, and j a is determined by the 
Gaussian part of the PDF. This gives a parabolic depen- 
dence on a. On the other hand, if a is large enough, 
the main contribution is due to a blob having a minimal 
possible volume which is of the order of L d . The decay 
exponent is fully determined by the probability to have 
such a blob and hence is a independent [SB]. Note that 
this picture implies that the exponential decay holds at 
t> lAdl^lnPe. 

If instead of a single blob one takes a spatially homo- 
geneous problem, this consideration should be slightly 
modified. At large t, initially uncorrelated blobs are 
brought close to each other because of the contraction 
along a certain direction, and then they overlap diffu- 
sively. Since the number of overlapping blobs is large, 
due to the central limit theorem it is rather d 2 which 
is inversely proportional to volume. Therefore, oc 
(exp (-a J2 Pi/2)). 

Formally, one should solve Eq. (1.2) with the initial 
condition i9(0,r) = i)o(r). The solution is 

/dk 
^-ptfo (W T (t)k) exp l-Q^K] , (3.2) 

Q(t) = K / dt'W^W-^t') [W{t)W- 1 {t')] T . (3.3) 
Jo 

From the qualitative arguments it is clear that the long- 
time asymptotic should be independent of the particular 
form of the distribution. We will take the simplest statis- 
tics, which is Gaussian with the pair correlation function 

(i?o(ri)0o(r 2 )>=x(ri2), X = Xo exp[-r 2 /(8L 2 )] . (3.4) 

This particular form is chosen for further convenience. 
In what follows we set L = 1. It is possible to generalize 
the calculation for arbitrary x and show that the results 
are independent of its form. 

To proceed, we introduce the generating function of 



%) = (apMM)]),, 



(3.5) 



Here we assume averaging over the statistics of a and 
the initial distribution of the scalar. The simplest part is 



to perform averaging o ver t he Gaussi an f ield i3q. To do 
this, we substitute Eq. (3.2) into Eq. (3.5) and using the 
expression for the characteristic function of a Gaussian 
random variable [tLTfl, obtain 



z(y) 



exp 



' 2 



dk 



; X (W T k)e 



Substituting %(&) = (87r) d / 2 xo exp(— 2k 2 ) and integrat- 
ing over k we obtain 



Z = ( exp 



Xo 



2 y/det I{t) 



(3.6) 



We used the fact that / = WW T + Q, which can be ver- 
ified by writ ing equation on WW T + Q and comparing 
it with Eq. ( |E4|) 

Using Eq. (|3.6|), one can find 



(m,0)\ a ) =C Q ((det/)- Q / 4 



(3.7) 



Here C a is a numerical constant. Equation (3.7) re- 
duces the problem to an averaging of powers of det / = 
Ilexp(2p l ), 



= C Q J d d p exp 



d 

a \ ^ 



(3.8) 



where V is the probability density function of p discussed 
above. In the large time limit this integral can be calcu- 
lated in the saddle-point approximation. The calculation 
is slightly different for d — 2 and 3. 



A. Two-dimensional case 



In d = 2 integral (3.8) should be calculated only over 
pi, since the distribution over p2 is stationary. The 
saddle-point equation is 

It is clear that p\ oc t. As long as one can use the 
quadratic expansion, which valid at least at small a, the 
solution of this equation is p\ = (Ai — aCu/2)t, hence 



a ( aCn 
7Q= 2 Al ~ — 



(3.9) 



At a > a cl - — — 2S' 2 (— \\) the val ue o f pi becomes much 
smaller than Ait; the integral ( |3.8| ) is determined by 
the boundary of the integration region, and therefore 
la = S^i— Ai), independent of a. 



The domain of validity of Eq. (3.9) depends on the 
value of the parameter Af. If it is much smaller than 
unity, we can use the quadratic approximation to S2 
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everywhere. This case effectively corresponds to the 
Kraichnan limit || . 

In the opposite limit Af > 1, the quadratic expansion 
of S 2 cannot be used for a > l/(fCn), and Eq. ( |3.9| ) is 
valid only at a <C l/(fC%\). The form of the intermedi- 
ate region is not universal, and depends on the particular 
form of S 2 and hence on details of the statistics of a. 



B. Three-dimensional case 

In d = 3 the result is similar to that of d = 2, though 
the consideration is slightly more complicated, due to the 
presence of an additional degree of freedom. There are 
two c ases to be considered. If A2 < 0, then V is given by 
Eq. ( 2.13 ) and the calculation is the same as for d = 2. 

If A2 > 0, both degrees of fre edom pi and p 2 are active 
and one should use the PDF ( 2.12 ). The saddle-point 
equations are 



dS 3 (x 1 ,x 2 ) 
dxi 



a 
2 



0, 



0S 3 (x 1 ,x 2 ) 
dx 2 



+ 77 = 0, (3.10) 



where x\ = (pi — \\t)/t and x 2 = (p 2 — X 2 t)/t. Again, 
the beginning of the curve is determined by the Gaussian 
part of V, and 7 Q is parabolic: 



la 



a 



:s:s 



(3.11) 



with Pl = (Ai + aC 13 /2)t and p 2 = (A 2 + aC 23 /2)t. 
At a larger than a cr calculated below, we have the a- 
independent behavior 



la = S 3 (-\i,~\ 2 ) . 



(3.12) 



Depending on the parameters, two different types of 
behavior can occur at a < a cr . First, it is possible that 
as a increases, p 2 will grow more slowly with t and at 
certain a will become much smaller than X 2 t. At larger 
a the integration over p 2 will be determined b y the re- 
gion p 2 <C X 2 t and one should replace system ( 3.10 ) by 
the single equation 



dS 3 (x 1 ,-X 2 ) 
dxi 



(3.13) 



In this case a cl - = -2 d 1 S 3 (x 1 , -X 2 )\ Xl= _ Xl ■ 

The other possibility is that at certain a the differ- 
ence pi — p 2 becomes much smaller than pi. 2. Because 



of the constraint pi > p 2 , for larger a integral (B.S) is 
determined by the boundary pi — p 2 of the domain of 
integration and the saddle-point equation becomes 



dS(xi,x 2 ) dS(xi : x 2 ) 



dxi 



dx 2 



0. 



(3.14) 



Then a cr = — [d±S(xi, x 2 ) + d 2 S{xi, x 2 )] at xi = — Xi 
and x 2 = — X 2 . Geometrically, the first case corresponds 



to elongated ellipsoids, and the second one to sheets hav- 
ing two largest dimensions of the same order. 

If the parameter Af is small enough, these changes of 
the regime occur within the Gaussian part of S 3 . Then 
one can perform a more detailed investigation. The first 
regime is realized if Ai > X 2 Ci 3 /C 23 and C 23 < 0. Then, 
at a > —2X 2 /C 23 Eq. ( |3.11 ) should be replaced by 



8 \Q 



r 2 



-C n 



-22 



Xi 



C12A; 
C: 



22 



\ 2 

A 2 



2Q 



22 



(3.15) 



If A2 = 0, Eq. ( 3.11 ) has no r egion of validity, and at 
a > one should use Eq. ( 3.15| ) which becomes 



a a (r cl * 



(3.16) 



In particular, within the Kraichnan model, Eq. ( 3.16] ) 
is valid for < a < a cr - Substituting Ai and Cy (see 
Appendix ^) one finds a cr = 4, 



3Da 



la 



(-!)- 



(3.17) 



for a < a cr and 7 Q = 3D for a > a ( 
different from the one obtained in Ref. 



Our result is 
, which coin- 
cides rather with Eq. (3.11). An exact solutio n for a = 2 
(see Appendix ^ and Ref. [||) supports Eq. ( 3.17 ). The 



reason for the discrepancy is the following. Despite the 
fact that p 2 tx ln(«;/D), it is impossible to ignore it com- 
pletely. If this were done, the anticorrelation between pi 
and p 2 , existing due to the incompressibility condition 
would lead to the growth of p 2 , thus making the calcula- 
tion inconsistent. 

The second regime takes place if C 23 > Ci 3 and 
C23 > C13A2/A1. T hen, starting from a — 2(Ai — 
X 2 )/(C 23 — C13) Eq. ( 3.11 ) should be replaced by an- 
other formula. Although the dependence on a is still 
parabolic, the coefficients are rather cumbersome, so we 
do not write this here. 



C. Gradients of the decaying scalar 

In the same manner one can consider the decay of the 
gradients of the scalar. In analogy, we can look for the 
correlation functions (|u;| a ), where a; = V$(i, r)\ r=(j . As 
in the case of single-point scalar statistics, these cor- 
relation functions decay exponentially in time. It was 
claimed in Ref. || that the decay law of the scalar and its 
gradient is the same within the Kraichnan model. Here 
we show that this is actually the case for arbitrary cor- 
related strain. Qualitatively it follows from the estimate 
that |V?9| ~ where I = exp(p^) is the smallest di- 
mension of the blob. As explained above, $ and I can 
be considered as independent, while the statistics of I is 
stationary. Thus the decay of the gradient is solely due 
to the change of the density of the scalar. More formally, 
one has 
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dk 

(2^V 



kcc'&o (W T {t)k) exp [-Qwkpkv] (3.18) 



Introducing the function Z(y) = (exp [i ( y, y )]) , and av- 
eraging it over the initial distribution (3.4) we obtain 
a formula similar to Eq. ( |3.6| ) Then, making a Fourier 
transform over y, we obtain the PDF of uj: 



( (det/) d / 4+1 / 2 exp 



VdetT 
Xo 



(w,Jw) 



(3.19) 



Considering this expression in the eigenbasis of the ma- 
trix /, we observe that (|w| Q ) ~ {\ud\ a ), since pd is 
smaller than the rest of the p is. Recalling that the dis- 
tribution over pd is stationary, we immediately obtain 
that 

(|V^(t,0)r)«((det/)- Q / 4 ) CT , 



which, due to Eq. (|3.7| ), gives the same law of decay. 



IV. FORCED TURBULENCE 

A. Single-point distribution of $ 

In this section we shall investigate the steady state dis- 
tribution of a passive scalar which occurs in the presence 
of a stationary source. For this purpose we introduce 
a random function (f>(t, r) on the right-hand side of Eq. 
(1.2), injecting blobs of the scalar with the character- 
istic size L. Due to the linearity of the problem, the 
scalar field at the moment t — is given by a superpo- 
sition of the scalar injected at earlier instants of time. 
Each realization of a can be characterized by a param- 
eter f* (cf. Ref. ||), such that the smallest dimension 
of blobs injected at t w — t* approaches r^a at t = 0. 
The ambiguity in the definition of t* is of the order of 
I Ad | _1 which is much smaller than the typical stretching 
time |A £ j|~ 1 lnPe. Considering the motion of the scalar 
injected at — 1„ < t < one may neglect diffusion, and 
the scalar is simply transfered along Lagrangian trajec- 
tories. On the other hand, as discussed in Sec. Ill, the 
contribution of the scalar injected at t < — t* is exponen- 
tially small. Thus t* separates diffusive and diffusionless 
regimes. One can write the following approximate for- 
mula: 



0(0, r)U 



dt(f>(t,0) 



(4.1) 



If the correlation time of the source is m uch smaller than 
t*, for a fixed realization of a integral ( |4.l| ) can be con- 
sidered as a Gaussian variable with zero average and the 
dispersion proportional to i„ , so that after averaging over 
(j>, for the single-point PDF H one obtains 



P(0) 



1 



exp 



t) 2 



(4.2) 



where Xo = J dt {cj>(t, 0)<p(0, 0)). The effective Gaussian- 
ity of the pumping has its limitations due to finite cor- 
relation ti me of 4> [analogously to the discussion of Eqs. 
( |2.7l) and ( |2.8| )1. Since we work in the comoving reference 
frame, this correlation time is very small, and hence only 
the tail will be affected. At the end of the section we 
discuss the implications of this. To proceed with formal 
calculations, we should specify statistics of <fc. Here we 
shall take <f> as a Gaussian field with the pair-correlation 
function 



where x( r ) is the same as in sec. III. The expression for 
the generating function Z = (exp(iyi?)) follows 



exp 



y 2 xo '•" 



dt' 



^oo\/ det/ ( ^') 



(4.3) 



Here I(t, t') is a matrix, satisfying Eq. Jl.4| ) with respect 
to t and the initial condition I(t',t') — 1. Naturally, the 
integration is performed over the initial time, summing 
up th e bl obs injected at different times. The integral in 
Eq. ( |4.3| ) gives the formal definition of i» entering Eq. 
(©: 



dt' 



-oo y/detl(0,t') 



(4.4) 



Intro duc ing the distribution function p(t*), we rewrite 
Eq. (|J) as 



dt* 



o V 27r Xoi* 



pit*] exp 



(4.5) 



Since t* is a functional of the whole trajectory p^(t, t'), 
one needs more information than contained in the simul- 
taneous distribution function of p. However, the follow- 
ing approximation, becoming exact at InPe — > oo, re- 
duces the problem to single-time statistics. We shall ne- 
glect the configurations for which the smallest dimension 
of the blob once reached r^if starts to grow. Then the 
realizations for which t* is larger than some T, and those 
for which the blob injected at — T has p 3 (0,~T) > r dif , 
are the same, leading us to the following formulas: 



p(t* 



d 



dpi / dp 2 V(U,pi,p 2 ) , (4.6) 

Ci* J-oo Jln( K /|A 2 |) 

Q POO POO P OO 

at * J-oo J-oo J1ii(k/|A 3 |) 



Here one should substitute PDF's (2.7) and ( |2.8| ), since 
t* is determined by the diffusionless regime. These equa- 
tions define nothing but the flux of the probability out of 
the region p 3 > ln(«;/|Ad|), once returns are disregarded. 
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Investigation of the above integrals shows that 
has the following properties. Its main body is concen- 
trated in the vicinity of t* = |Ad| _1 lnPe, and has a 
width of the order of VhiPe. On the other hand, its 
tail t* |Ad| _1 InPe decays exponentially, 



p(t*) oc exp(— ci*) , 



(4.7) 



where c is equal to S , 2(— Ai) in d = 2 and Sa(— Ai,— A3) 
in d = 3. The intermediate region is not universal and 
depends on the details of a. Note that since Eq. (4.7) 
gives the probability that at large time a blob has not 
yet decayed, therefore c is equal to the limiting value of 
7„ (see sec 



III) 



This information allows one to calculate the probabil- 
ity distribution function V{~d). If # <C Ir±(|Ad//«), it is the 
central peak of p(t*) that determines the scalar PDF: 



V = 



2ttxo hi P° 



1/2 



exp 



|A d |^ 2 
2 Xo lnPe 



(4.8) 



As we increase 1?, at i9 > InPe the details of the dis- 
tribution of t* become important. The Gaussian regime 
[Eq. ( |4.gD ] will turn into some nonuniversal asymptotic. 
Nonetheless, at $ ^> InPe, due to Eq. (4.7) the univer- 
sality is restored: 



V tx exp 




(4.9) 



For a 5-correlated a one can find the complete func- 
tion V{°0). Then $2,3 is Gaussian (Appendix |A]) and the 
result can be found in the saddle- point approximation. 
In d = 2 the result coincides with that of Rcfe. P , |il||ll| . 
In d = 3 we obtain the formula (cf. Refs. [^o| , ^T]| ) 



\nV oc -3 



/Dd 2 
In 2 PeH InPe 



2xc 



for 1 1? | < 4^/xo/DlnPe, and 



ln-p oc 61nPe-4W31n ii Pe 



3Dd 2 



otherwise. The change of the regime is related to the 
fact that the two dimensions of the contributing blobs 
start to be equal to Td\£- This result is different from the 
one presented in Ref . [[l0| . The difference can be qualita- 
tively explained as follows. In our case the structures of 
the scalar making the main contribution to the PDF are 
columns, with the two smallest dimensions of the same 
order. They appear because of the anticorrelation origi- 
nating from the incompressibility condition: for t* larger 
than a mean value, P3 should decrease slower than X^t, 
which by virtue of the anticorrelation leads to a decrease 
of p2 faster than A2^. Staring from a certain value of 
1?, both pi and p$ decrease at the same rat e (an analo- 
gous phenomenon is described in subsection [II B). This 



structure is different from the ansatz proposed in Ref. 

Equation (4.2) should be modified if <f> has a finite cor- 
relation time [see the discussion leading to Eqs. (2.7] 
and (O) ). That is 



T 3 — ( —^== ex P 



where f(x) deviates from x 2 /(2xo) at s > V T 0- This 
may affect only the tail of V($). If the parameter t^^/xqc 
is much smaller than unity the tail is determined com- 
pletely by the region where / a i 2 , and one obtains 



the asymptotic result (4.9). Conversely, at 



the form of / should be accounted. Nevertheless, one 
can easily check that the exponential tail survives with a 
decrement depending on the form of /. 



B. Gradients 



Here we briefly consider the statistics of the scalar gra- 
dients. Within the Kraichnan model the problem was 
solved in Ref. |jj . From the qualitative picture presented 
there, one can conclude that the PDF is determined by 
the short-time fluctuations of er, and hence is nonuniver- 
sal. The following considerations support the con cl usion . 

In a wa y, sim ilar to the one leading to Eqs. (3.6), ( 3.1£ , 
(4.2, and (4.3) one can find 



Z{y) = ( exp 



, y/detl(0,f) 



(4.10) 



Analogously to the case of the scalar density, the gradient 
field at t = is given by a superposition of contributions 
of blobs injected at earlier moments of time. We observe 
that the contribution of each blob into ui 2 = (V$) 2 is de- 
termined by two factors: the value of the scalar density 
Xo(det I) -1 / 2 and the inverse size of the blob contained in 



I a p. Not all the blobs make a contribution to Eq. (4.10) 



Indeed, the size of the blob s inje cted at |f'| <C |t*|, (where 
f * is defined in subsection IV A ) is much larger than rdif 
and, therefore the value of the gradient will be small. On 
the other hand, the scalar injected at \t'\ 3> |t*| has an ex- 
ponentially small density and hence does not contribute. 
Thus the distribution of lu is determined by the blobs in- 
jected at t' ~ — t* which have the minimum possible size 
provided the diffusion is still ineffective at t = 0. 

Each realization of a can thus be roughly characterized 
by two relevant parameters. The first one is the lateral 
dimension I of the thinnest blobs, for which the diffusion 
can still be neglected at t < 0. This is related to the 
very last stage of the evolution, when blobs of the small- 
est size of the order of rdif may undergo a strong rapid 
contraction, increasing the gradient without dissipating 
the scalar. Let us stress that the fluctuation should be 
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short lived in order to suppress the diffusive spreading of 
the particles. 

The other parameter is the duration Tq of the injection 
stage for these blobs, showing how many blobs approach 
rdif at t f» 0. There is no average strain during this pe- 
riod, so that blobs injected at — f* — To < t' < — 1» all 
have a size of the order L at t s» — i*. Since at t > — f* 
the blobs move in the same velocity field, they all have 
approximately the same size at t = 0. Formally, the num- 
ber of relevant blobs is expressed by the formula [c.f. Eq. 
©] 



-t. 



dt' (/>(?, 0). 



— Tq— t* 

Writing the estimation for the gradient w 
replace Eq. ( |4.10| ) by Z = (exp(%??/O) 0) , 
over </>, we obtain 



(4.11) 

#// we can 
Averaging 



exp 



y 2 xoT 

2l 2 



To, l 



Since the injection stage occurs at \t\ > |t»| w 
|A<j| -1 lnPe 3> t, the fluctuations determining To and 
I are independent |7|||| . 

On the average To ~ | | 1 and I ~ rdif, so that 
(uj 2 ) ~ Xo/(A^" 1 rj if ). Nevertheless, studying the tail of 
the gradients PDF, it is necessary to take into account 
the large deviations of these parameters. The probability 
of a large value of To is relate d to configurations of small 
strain (see Sees. QI and IV A ), and decays as exp(— cT ). 
Writing 



V(u) ~ (expl-LJ 2 l 2 /(2 Xo T )}) 
one can average over Tq and find 



Tod 



exp[ 



\lo\1(2c/xo) 1/2 \ 



To average over I one notes that the tail I <C r<jif of the 
probability distribut ion functio n of I is related to the tail 
ofVst(pd) [see Eqs. ( 2.11 - 2.13 )] via I = exp(pd)- Indeed, 
both are determined by the probability of a strong and 
rapid contraction from rdif to I <C rdif. Hence 

JdlV st {\nl)exp [-\lu\1{2c/xo) 1/2 ] (4.12) 

Within the Kraichnan model lnV s t oc —l~ 2 (see Ap- 
pendix A) and the result lnV(u) oc -|cj| 2 / 3 of Ref. § 
easily follows. In general the fluctuations of the smallest 
dimension take place at times of order r near t = 0, and 



therefore are related to the single-time distribution of a, 
which is nonuniversal. 

We conclude that the gradient statistics is nonuniver- 
sal, and cannot be predicted unless some specific informa- 
tion is supplied. For example, if the distribution of I falls 
off very fast at small the distribution of uj will have 
an exponential tail at very large to. In particular, this 
can explain the results of numerical simulations P,^9[, 
where a cutoff can be due either to the grid step or im- 
posed by hand || . If the tail of V s t behaves according 
to ln[7 7 s t(2)] oc — l~ a , the tail of V(u>) has the stretched 
exponential form: ln[P(a>)] oc — |cj| q /( q+1 ). 



V. CONCLUSION 

We considered a passive scalar advected by a random 
large-scale velocity field. Our purpose was to establish 
the degree of universality of the scalar statistics for ar- 
bitrary an correlated velocity. The investigation can be 
reduced to the statistics of different Lagrangian charac- 
teristics of the smooth flow. In the limit of large a Peclet 
number, part of the relevant information is contained in 
the long-time asymptotic properties of the Lagrangian 
statistics, which is shown to possess a universal form. 
The scalar quantities related to long-time evolution thus 
manifest universal statistical features. Generally, these 
are the central part and the tail of the corresponding 
PDF. We considered several particular examples of such 
quantities: the decay of the scalar density and its gradi- 
ent, and the scalar density in the forced case. Conversely, 
the statistics of the gradients in the forced case requires 
information about short-time fluctuations of the velocity, 
and is thus sensitive to its details. 

The application of the Lagrangian statistics estab- 
lished here is not only restricted to the above examples. 
One can slightly modify the procedure to consider many- 
point correlation functions of the scalar, say the PDF of 
the scalar difference at two points flCfl, correlation func- 
tions out of the convective interval |12|, and other prob- 
lems. A similar scheme could be applied to other passive 
quantities, like vectors Q and tensors 
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APPENDIX A: THE KRAICHNAN MODEL 



Gaussian with the pair-correlation function jjl 



Here we consider Eqs. (2.2,2.3) within the Kraichnan 
model, when matrix a has zero correlation time and is 



{a a ff(t)a,j, u (0))=D[(d + l)SaftS^ v — S cl pS llv — 5 av 5p l ^S{t) . 
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The tensor structure is fixed by incompressibility condi- 
tion. Zero correlation time allows to write the Fokker- 
Planck equation for the probability distribution of R and 
A. Integrating out the angular degrees of freedom, one 
can see that the equation obtained is equivalent to the 
following Langevin dynamics || 



dtPi 



Dd 
~2~ 



d 

E 



coth(2pi - 2pj) + Ci + - exp(-2pi) 



where ^ are random Gaussian delta-correlated processes 
with the following correlation functions 

= <V(* ~ . Cu = D(d5 lJ - 1) . 

Let us now consider a typical evolution of the eigenvalues. 
At t = all the eigenvalues pi are equal to zero. Then, 



P(t,pi,p 2 ) 



1 



exp 



(Pi - Dtf 
2Dt 



0(pi)S{pi + P2 ) 



during short initial period of time, all p start to differ. 
We can always arrange p so that pi > p 2 > .. > p d - We 
observe that then the ballistic terms coth(2pi — 2pj) 
are arranged in the same order so that the eigenvalues 
will continue to separate and at t 3> D^ 1 the following 
inequalities will hold pi ^> p 2 ^> ■ • ■ 3> pd- If this is the 
case, we can substitute the hyperbolic cotangents by ±1 
and get the following equations 



dtPi = Xi + & 



- exp(-2p 4 



), Xi 



Bd (A 



2i + l). 



This simplified dynamics can be easil 
probability density functions of p |l3| 



turned into the 



(Al) 



'Pit, Pl,P2,Pz) 



2V3irDt 



exp 



[i Pl - 3Dt) 2 + (p! - 3Dt) P2 + pi] 



3Dt 



Oipi - P2)0ip 2 - p 3 )S( Pl +p 2 + pa) . (A2) 



As explained in the main text these formulae arc valid 
at times 1/D -C t < 1/ 'D\n{DL 2 / '«) and the form of 
the PDF near the boundary can be modeled by the step 
function. At times t 3> l/-Dln(Z?L 2 /K) one should also 
calculate the stationary PDF of pd- They are readily 
found from the one-dimensional Fokker-Planck equation 



P 8t (p 2 ) = ^exp(-2p 2 -^e-^ 



PstiP3 



1 / K \ 3 / 2 / K 



(A3) 
(A4) 



APPENDIX B: SMALL r-EXPANSION 

In this section we assume that the correlation time 
t of (j is small, namely Dt -C 1, where L> = 
(tr - / dtu 7 (O)u(t)) characterizes the amplitude of the fluc- 
tuations of cr. We investigate the effect of finite correla- 
tion time on the Lyapunov spectrum. In d = 2 there 
are no essential changes with respect to the ^-correlated 
case since both Ai and X 2 = — Ai get small corrections 
in t leaving all qualitative features unchanged. However 
in d = 3, one can ask whether A2 will shift from its zero 
value at r = 0, and whether the correction is positive or 
negative at finite r. We demonstrate that already the 
first order correction in r leads to a generally non-zero 
value of A2 in d = 3, which can be both positive and 
negative. 

Here it is more convenient to parameterize the an- 
gular degrees of freedom of WW T by the eigenvectors 
et instead of matrix R (see Eqs. (2.2, |2.3| )), given by 
ef = Ri a . If the eigenvalues are separated, that is 



Pi 3> p 2 3> P3 , the equations for pi and corresponding to 
it eigenvector ei decouple 

d t pi = (ei,crei) , 9 t e 1 = ae 1 -e 1 (e 1) ae 1 ). (Bl) 

The same is true for p% and 

dtP3 = (e 3 ,cre 3 ) , d t e 3 = -a T e 3 + e 3 (e 3 ,ere 3 ) . (B2) 

This system implies that under the transformation a — » 
— a T the eigenvalues are transformed as Ai,3 — » — A34. 
In the calculation it is convenient to deal with symmet- 
ric matrices, which is achieved by decomposing a into 
symmetric and antisymmetric parts a — s + uj and intro- 
ducing ei = il/n with 



so that 



d t M = ojM , M (0) = 1 



Ai = ((n, In)) , 9 ( n = sn — n(n, sn) 



with s = M T sM . To find the first order correction to 
Ai we integrate the above differential equation from to 
t, and then iterate the obtained expression once. After 
averaging over the directions of n(0) we find 



Ai = - / dh 
5 Jo 



tr{s{t)S(ti)) + - / dt 2 tr (s{t)s(ti)s(t 2 )} 
' Jo 



where t 3> t. The expression for — A 3 is obtained by 
changing s — > — s in this formula. Then, using X 2 = 
— Ai — A 3 , we find 



A - 12 



Jo 



dtidt 2 tr(s(t)s{ti)s(t 2 )) 



which is generally non-zero and has no definite sign. 
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APPENDIX C: PAIR-CORRELATION 
FUNCTION 

In this appendix we calculate the time decay of the 
pair correlation function f(t,r) = (&(t, r)$(t, OY) within 
the Kraichnan model. It satisfies the equation 0] 

d t f = JC aP {r)V a Vpf + 2kV 2 / , f\ t=0 = X (r) , 



K 



a/3 



D 



Making the Fourier transform over r and passing to 
spherical coordinates, we get 



d T f = k 2 d 2 f+(d+l)kd k f-ek 2 f. 

Did -I) 4k 
-t , e = 



(CI) 



2 ~ D(d-l) 

Next, making the Laplace transform over r we obtain 

k 2 f" + (d+ l)kf -(E + ek 2 )f = -x(fc) • (C2) 

Two linearly independent solutions of the homogeneous 
equation are expressed via the modified Bessel functions 
of the order v = y/E + gP/4. The branch of the square 
root should be picked up so that it has a cut along 
the semiaxis E < —d 2 /A and takes positive values at 
E > — gP/4. Using these func tion s one can find the Green 
function g(E,k,k') of Eq. (C2) satisfying the correct 
boundary conditions 

g(E, k, k') = k-^k'^ 2 - 1 \l v (k^ 2 )K v (^ 2 )0^ - k) 



+K v {ke 1/2 )I v {k' 'e 1/2 )9(k - k') 



(C3) 



with the solution of Eq. ( |Cl| ) given by f(t, k) — 
fdk' g(t, k,k')x(k')- Next we should make the inverse 
Laplace transform of g(E, k, k') 



-i />b-\-ioo 

g(r,k,k') = — dEe ET g{E,k 1 k') 
tni Jb-ioo 



(C4) 



Here b > is arbitra ry. One may deform the integra- 
tion contour in Eq. (C4) until a singularity of the in- 
tegrand is encountered. The first singularity appears at 
E = —d 2 /4, which is the branch point of v. Therefore, 
the integration should be performed along the real axis 
at — oo < E < — d 2 /A on both sides of the cut. Making 



the change of variable E = — . 

g(r,k,k') = 2_ k -d/2 yd/2-i e - 



d /4 we obtain 



d t/4 



< I dxe 
'a 



r xsmh(irx)K lx {ke 1/2 )K tx (k'e 1/2 ) 



If one is interested in the single-point statistics, one 
should integrate this expression over k 



dkk d ~ l g{T,k,k') 



od/2-l -d/i , „ 
z e yd/2-l n -d 2 r 



poo 

x / dxe~ x2T xsmh(-Kx)K ix (fc'e 1 / 2 ) 



r/4 



d ix 
4 + Y 



Al large r the integral is determined by a narrow vicinity 
of x — 0. After a simple calculation we obtain 



($ 2 (t)) 



CPe' 



d/2 



(Dtf/ 2 



exp 



d 2 (d-l)Dt 



(C5) 



where C is a x-dependent constant. 
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